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Abstract 

We investigate the dynamics of the spin-boson model when the spectral density of 
the boson bath shows a resonance at a characteristic frequency but behaves Ohmi- 
cally at small frequencies. The time evolution of an initial state is determined by 
making use of the mapping onto a system composed of a quantum mechanical two- 
O ■ state system (TSS) which is coupled to a harmonic oscillator (HO) with frequency 

O. The HO itself is coupled to an Ohmic environment. The dynamics is calculated 
by employing the numerically exact quasiadiabatic path-integral propagator tech- 



nique. We find significant new properties compared to the Ohmic spin-boson model. 



By reducing the TSS-HO system in the dressed states picture to a three-level system 
for the special case at resonance, we calculate the dephasing rates for the TSS an- 
alytically. Finally, we apply our model to experimentally realized superconducting 
flux qubits coupled to an underdamped dc-SQUID detector. 
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1 Introduction 



The spin-boson model is one of the central paradigms in theoretical physics 
[1,2]. The quantum mechanical two-state system (spin 1/2 or in today's mod- 
ern language qubit) interacting with a bosonic bath, is the simplest possible 
model to describe the effect of an environment on constructive and destructive 
quantum interference. It allows in a rather comprehensive way to investigate 
decoherence and damping on the quantum system imposed by a phenomeno- 
logical environment. Hence, it has found numerous applications ranging from 
electron transfer to quantum information processing. 

The characteristics of the harmonic bath in contact with the two-state sys- 
tem is captured by the spectral density J(uo) of the bath. A well studied 
special case is the spin-boson model with an Ohmic spectral density, where 
J{uj) oc uj (up to some cut-off frequency). This implies that the quantum sys- 
tem is damped equally at all frequencies, which is the case for many physical 
systems. A prominent example in the field of condensed matter physics are 
single-charge effects in micro- and nanostructured systems in which the electro- 
magnetic environment can be modeled by an impedance of pure Ohmic form. 
The more general case of J{w) having a power-law distribution, i.e., J{oo) oc oo s 
(0 < s < 1 sub-Ohmic and s > 1 super-Ohmic bath) also finds wide-spread 
applications, for instance when the TSS couples to an electromagnetic envi- 
ronment represented by a RC-transmission line. For these specifications, the 
time constants with which a quantum superposition of the TSS decays can be 
calculated in analytic form by using real-time path integral techniques [1,2]. 
Another limit to be mentioned in this context is the case of a single environ- 
mental mode at a specific frequency Q, i.e., J(ou) oc 5(u — Q), which is for 
instance realized by an electromagnetic environment with a pure inductive 
impedance. 

In this work, we consider the spin-boson model for the case of the bath spec- 
tral density having a distinct resonance at a characteristic frequency Q of the 
order of the TSS level splitting, and behaving Ohmically a low frequencies 
(peaked spectral density with Ohmic background). This model currently re- 
ceives growing interest in the context of quantum computing with condensed 
matter systems [3], especially for superconducting flux qubit devices [4,5,6,7]. 
In these devices, the state of the TSS, i.e., the direction of the magnetic flux, 
is read out by a dc-SQUID which couples inductively to the qubit supercon- 
ducting loop. At bias currents well below the critical current I C) the SQUID 
can be modeled as an inductor Lj. An additional shunt capacitance C s creates 
an on-chip environment which improves the resolution of the dc-SQUID. The 
resulting impedance of this dc-SQUID environment can display a pronounced 
resonance at the characteristic frequency of the LC-circuit (see below for de- 
tails). 
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As an alternative but equivalent point of view, a TSS coupled to a harmonic 
bath with a peaked spectral density represents an effective description of a TSS 
coupled to a damped harmonic oscillator of characteristic frequency Q. In this 
perspective, this model has recently been used to describe solid state imple- 
mentations of a qubit coupled to a resonator. Various schemes of measurement, 
entanglement generation, quantum information transfer and interaction with 
a nonclassical state of the electromagnetic field have been proposed [8,9] and 
dephasing effects have been addressed [9]. Moreover these devices, besides be- 
ing promising candidates for realizing qubits in a possibly scalable quantum 
computer, are also very interesting objects to display fundamental quantum 
physical properties, like macroscopic quantum coherence [10], Rabi oscillations 
[6,11,12], Ramsey interference [6,13], etc. 

The spin-boson model with a peaked spectral density of the bath has been in- 
vestigated theoretically in Refs. [14,15]. In Ref. [14], it has been demonstrated 
that an external driving field can be applied to slow down decoherence by 
moving the TSS out of resonance with the HO due to dressing of the TSS by 
the driving field. Recently, Kleff et al. [15] have calculated for the static case 
at zero temperature dephasing times by applying the numerical flow equation 
renormalization method (see below for further remarks). 

To explore the more realistic finite-temperature dynamics of the spin-boson 
model with a resonance peak at frequency Q in the spectral density of the 
bath, we make use of the already mentioned exact mapping to a TSS+HO 
model [16]. By this, we obtain a model consisting of the TSS which couples 
to a harmonic oscillator with characteristic frequency Q, which itself couples 
to a harmonic bath with Ohmic spectral density. Throughout this work, we 
concentrate on the case of a weak coupling of the oscillator to the Ohmic bath, 
but the TSS-HO coupling is kept arbitrary. This choice of the parameter region 
is due to the experimentally given situation of dc-SQUIDs which are typically 
underdamped. This mapping allows to apply the numerically exact method of 
the quasi-adiabatic propagator path-integral (QUAPI) [17] for the TSS+HO 
being the the central quantum system. Ay additionally tracing out the HO 
degrees of freedom, we determine the time-evolution of the reduced density 
matrix of the TSS. We extract the dephasing and relaxation time constants for 
the decay of a localized state of the TSS. In order to elaborate the difference, we 
compare our numerical results to the well-known analytical expressions for the 
rates of a TSS coupled to a structured harmonic reservoir evaluated to lowest 
order in the TSS-reservoir coupling [18]. As it turns out, these weak-coupling 
rates deviate considerably (up to a factor of 50) from the numerically exact 
decay rates for level splittings around the oscillator frequency. This shows that 
if the frequencies of the TSS and the HO are comparable, the decay rates are 
determined by the full frequency spectrum of the peaked environment, and 
not only by its low frequency weak Ohmic backround. Next, we concentrate 
on the specific case of the TSS being exactly in resonance with the oscillator 
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and establish a simple three-level description of the TSS-HO system. Since 
we are interested in weak system-bath coupling, we formulate the Redfield 
equations for this three-level system (3LS) which can be solved analytically. 
The resulting dephasing rates agree well with the numerically exact values 
of QUAPI. Finally, we apply our general model to experimentally realized 
superconducting flux qubit devices. 



2 The model and the mapping to a TSS+HO Hamiltonian 

To setup the model, we consider the Hamiltonian of a two-state system being 
coupled to a bath of (non-interacting) harmonic oscillators, i.e., 




Here, b k and b\ are the annihilation and creation operators of the k— th 
bath mode with frequency uJk- The spectral density of the bath modes is 
assumed to have a Lorentzian peak at the characteristic frequency Q, but be- 
haves Ohmically at low frequencies with the dimensionless coupling strength 
a = lim^o J{uj)/2u. The Lorentzian-shaped resonance has a width which we 
denote as 7 = 2itkVL. To be specific, we assume a spectral density of the form 

J(u) = E \l5(u - Q k ) = (fl2 _ J™+ {2kkuQ)2 -2a,. (2) 

For our purposes, it is convenient to aquire a different viewpoint as it will 
become clear below. Following Ref. [16,14], there exists an exact one-to-one 
mapping of this Hamiltonian onto that one of a TSS coupled to a single HO 
mode with frequency fl which itself interacts with a set of (non-interacting) 
harmonic oscillators having an Ohmic spectral density with the dimensionless 
damping strength k. Thereby, the interaction strength between the TSS and 
the HO is given by g. The corresponding Hamiltonian reads 

5 = -^ - ^-a z + hga z (B ] + B) + hVtB ] B 

+(£t + B) ]T hu k (b\ + b k ) + Y, Tuv k blh . (3) 
k k 

Here, B and B^ are the annihilation and creation operators of the localized 
HO mode while the b k and the b\ are the corresponding bath mode operators. 
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The spectral density of the continuous bath modes now becomes Ohmic, i.e., 
^Ohm(^) = u l s ( u ~ u k) = Kue - "^" , (4) 



where we have introduced the usual high-frequency cut-off at u c (We note that 
throughout this work, we use the value u c = 10A ). The relation between g 
and a follows as 

g = ^- « a = 8 K ^. (5) 



Since for our cases of interest, the damping of the HO is small (re <C 1), we 
will denote the TSS+HO-system as the central quantum system which is itself 
weakly damped. However, the coupling between the TSS and the HO can still 
become large. In the following, we will evaluate the reduced density opera- 
tor pTss+Hoif) = tr h&th e iHt / h W{Q)e- iHt / h for the TSS+HO-system, where the 
bath degrees of freedom have been traced out. W(0) denotes the full density 
operator at initial time t — 0. Throughout this work, we assume a factorizing 
initial preparation where all three parts of the total system are decoupled and 
the coupling is switched on instantaneously at t = + . In particular, we choose 
the bath being at thermal equilibrium at inverse temperature (5 = 1/ksT as 
well as the localized mode being thermally distributed. This implies 

e -f3H H o p-@H bath 

W(0) = p TSS (0) ® — ® — , (6) 

^HO ^bath 



where Z denotes the partition function of the corresponding subsystem. The 
further reduction to the reduced density operator p(t) of the TSS alone is easily 
done by tracing out the HO degrees of freedom, i.e., p(t) = tr HO prss+Hoif)- 



3 Observables 



Having calculated the reduced density matrix p(t) of the TSS, we can di- 
rectly calculate the observables of interest. In this context, we are interested 
in the dephasing and relaxation rates of the TSS dynamics. Hence, we prepare 
the TSS in a symmetric superposition of its energy eigenstates. The relevant 
dynamical quantity then is the symmetrized correlation function 

S,(t) = \(a z {t)°M + ^(0K(0) - {*,)l t+ , (7) 
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where cr z (t) = e iHt / h (x z e~ tm / h and where the equilibrium population (o" 2 )^ c + 
(the subscript + indicates that the dynamics has been calculated with positive 
bias e > 0) has been subtracted. By assuming a factorizing initial preparation, 
the correlation function can be expressed in terms of the population difference 
P+(t) = {<y z )t = P s (t) + P a (t). It follows 

S z (t) = P s (t) + Poo t+ [P a (t) - Poo !+ ] , (8) 



where P s / a (t) are the symmetric/antisymmetric parts of P+(t) with respect to 
the sign of the bias and Poo,+ = +• This quantity vanishes at long times 

which ensures that the Fourier transform exists. By calculating the population 
P-(t) with the sign of the bias inverted (e < 0), we obtain P s / a (t) according 
to 



P.(t) = [P + (t) + P4t)]/2, 

P a (t) = [P + (t)-P_(t)}/2. (9) 

For the unbiased case e = 0, we have S z (t) = P s (t) = P(t). Moreover, the 
initial preparation prssi®) of the TSS in an equally weighted superposition of 
ground and excited state implies P+{0) = P_(0) = 1. Since S z (t) vanishes at 
t — > oo, it readily can be Fourier transformed, i.e., 

oo 

S z (lu) = 2 J dt cos LutS z (t) . (10) 



To extract the decay rates Tj, the frequencies uoi and the spectral weights a* 
associated to S z (u), we approximate S z (t) as a sum of exponentially decaying 
sinusoids, i.e., 

s z(t) = a i e ~ Vlt COSCUit , (11) 



where J2i a i — 1 • It h &s the Fourier transform 



S z (u>) 



Tj + (ui + u) 



+ 



(12) 



By a standard fitting procedure, we fit the exact numerical results to this 
function to extract the decay rates Tj (HWHM of the Lorentzian peak), the 
frequencies oj,i (position of the peak) and the spectral weights a*. 



6 




Fig. 1. Example of the dynamics for the symmetric case e = 0, where the oscillator 
frequency is in resonance with the TSS frequency, i.e., O = Ao- Parameters are: 
g = 0.18A , k = 0.014(-» a = 0.004), k B T = 0.1£A . QUAPI parameters are 
M = 12, K = 1, At = 0.06/A . 

4 The numerical ab-initio technique QUAPI 

In order to calculate the reduced density operator, we use the iterative ten- 
sor multiplication scheme derived for the so-called quasiadiabatic propagator 
path integral (QUAPI). This numerical algorithm was developed by Makri 
and Makarov [17] within the context of chemical physics. Since its first ap- 
plications it has been successfully tested and adopted to various problems of 
open quantum systems, with and without external driving [17,19,20,21]. The 
details of this technique have been discussed to a great extent in the litera- 
ture [17,19,20,21] and we only mention briefly how to adopt it to our specific 
problem. 

We use QUAPI to calculate the reduced density matrix p(t) and iterate until 
the stationary state of P±(t) has been reached. For the biased case, we run 
the iteration with both positive and negative bias to calculate P±(t) (e > 
and e < 0). Having transformed to P s / a (t), see Eq. (9), we determine S z (t) in 
Eq. (8) and Fourier transform the result to determine S z (cu) according to Eqs. 
(10,). Then, a fit to Eq. (12) determines the dephasing and relaxation rates 1^ 
(HWHM of the Lorentzian peak), the frequencies Ui (position of the peak) and 
the spectral weights a*. A typical result for the symmetric TSS (e = 0) being 
in resonance with the HO (A = fi) is shown in Fig. 1, see also Section 9.2 
for a more detailed discussion of this specific case. The inset of Fig. 1 depicts 
the time evolution of P(t). The Fourier transform displays two characteristic 
peaks at uj = uj\ and uj = uj 2 . The width of these peaks can be read of the 
fitted Lorentzians, see dashed lines. 

QUAPI uses three free parameters which have to be adjusted according to 
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Fig. 2. Main: Fourier transform S z (lo) for two different memory lengths 
K = 1 and K = 2, see text. Parameters: e = 0,g = 0.18Ao, 
tt = 0.75A , k = 0.014(-> a = 0.006), k B T = 0.1M . Remaining QUAPI parame- 
ters are M = 12, At = 0.06/Ao- Inset: Dephasing rate T 2 for different QUAPI-time 
steps At for two different values of Q. The position of the maxima determine the 
choice of At for QUAPI, see text. 

the particular situation. Typical physically reasonable situations allow for a 
tractable choice of parameters. 

(i) The first parameter is the dimension M of the Hilbert space of the quan- 
tum system which is in our case the composed TSS+HO-system. With the 
dimension of the TSS being fixed, we truncate the in principal infinite dimen- 
sional oscillator basis to the six lowest energy eigenstates leading to M — 12. 
We have checked convergence also with M — 16 and M = 20 and found the 
results unchanged. 

(ii) The second parameter affects the length of the memory which is induced 
by the environmental correlations. This memory time has the length of K 
times the time step of the iteration. Here, we deal with weak Ohmic damping 
(k <C 1) of the HO and temperature in the experimental devices is typi- 
cally larger than damping. (This in principle would even allow for a Marko- 
vian approximation.) In our case, this permits to choose K small. We fix 
K — 1. We have also ensured convergence by choosing K = 2, but the re- 
sults changed only by less than 10 %. This is depicted in Fig. 2. The results 
for peak 1 are for K = 1 : I\ = 0.0177, wi = 0.730, ai = 0.197 (in units 
of A ) while for K = 2 : T t = 0.0166, u x = 0.718, a t = 0.181. For peak 2, 
we find the values K = 1 : T 2 = 0.0310, u 2 = 1.074, a 2 = 0.913, while for 
K = 2 : T 2 = 0.0282, lu 2 = 1.073, a 2 = 0.917. The sum rule a x + a 2 = 1 is 
fulfilled with an accuracy of approximately 10%. 

At this point, it becomes clear why it is essential for QUAPI to use the 
mapping to the TSS+HO model. For the dynamics of the spin-boson model 
with the peaked spectral density in Eq. (2), the autocorrelation function 
L(t) = (£(£)£(0)),£ > 0, of the fluctuating force £(£) of the bath [1] decays 
slowly with time in an oscillatory way. This can be readily understood since 
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L(t) is qualitatively the Fourier transform of J(u). In other words, the peaked 
environment possesses a longer correlation time compared to the Ohmic bath. 
This is also the reason why a standard weak-coupling approach fails to de- 
termine the dephasing rates (see Sec. 9.1 for a detailed comparison between 
both approaches). Since the long correlation time of the peaked bath would 
require a long memory time KAt with K ^> 1, the direct treatment of the 
TSS coupled to the peaked bath with QUAPI is still beyond the today's com- 
putational resources. Only the corresponding Ohmic autocorrelation function 
decays fast enough with time to be truncated after a relatively short memory 
time (at any finite temperature T). It is computationally more favorable to 
iterate the HO dynamics as a part of the quantum system and not as a part 
of the bath. 

We mention that for the comparison with the three-level approximation (see 
Section 9.2 below), which is expected to be good at low temperatures, we used 
the parameters M = 6 and K = 3, since in this case the bath-induced corre- 
lations have a weaker decay and hence require a longer memory time. Also, 
the lower temperature allows to restrict to less basis states, 
(iii) The third QUAPI-parameter is the time-step At. It is fixed according to 
the principle of minimal sensitivity [20] of the result on the variation of At. 
The results which permits to decide on At are shown in the inset of Fig. (2). 
The position of the maxima determine the choice of At. In our case, we fix 
At = 0.06/Ao- As one can observe from Fig. (2), the position of the maximum 
is not sensitive on the change of the parameter Q. 



5 Varying the oscillator frequency f2 

In this section, we investigate the effect of varying the frequency Q of the 
resonance of the spectral density J(oo) in Eq. (2), or, in other words, the 
frequency of the harmonic oscillator. For simplicity, we restrict here to the 
symmetric TSS e — 0. For the biased TSS, the qualitative behavior is similar 
if one replaces A by the effective level splitting A b = \J Aq + e 2 for the biased 
case. 

5. 1 Exact resonance Q = A 

The dynamics of the TSS being in resonance with the HO is depicted in Fig. 1 
for a typical set of parameters (see caption). One observes a decay of P(t) with 
two frequencies and two damping constants. This can be readily understood 
by considering the undamped TSS+HO-system, i.e., k — 0. To diagonalize 
the corresponding Hamiltonian, we can use the rotating-wave approximation 
(RWA), which is appropriate at resonance. Treating the interaction g in first 
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Fig. 3. S z (oj) for two values of the oscillator frequency Q. Parameters are: 
£ = 0, g = 0.07A , k = 0.014, k B T = 0.1hA . 



order perturbation theory, we obtain the spectrum as a ladder of energy-levels. 
The groundstate \g,0) (g denotes the groundstate of the TSS while denotes 
the index of the HO mode) is separated by nhVt from the higher lying n-th 
pair of excited states which are almost degenerate, but split by HA n = Ihg^fn. 
At rather low temperature, only the three lowest energy eigenstates, i.e., the 
ground state \g,0) and the first (n = 1) pair of excited states are relevant. 
This yields to two peaks in the spectrum at u 1 = A — g and at u 2 = A + g 
with the distance A 1 = 2g which appears in Fig. 1. The two peaks have almost 
the same height, i.e., the same spectral weight cq ps a 2 - 



5.2 Negative detuning Q < A 



If the HO frequency is detuned from the TSS level splitting such that Q < A , 
the dynamics of P(t) also contains two frequencies but with different spectral 
weights. This is shown in Fig. 3 for Q = 0.75Ao- Also, a perturbative analysis 
of second order in g (the first order vanishes exactly), readily explains the 
features: again restricting to the three lowest energy eigenstates, we find that 
the first excited state \g, 1) is separated from the ground state \g, 0) by lo\ = 
^ — a^-q. 2 ~ ^ — Ao-n > wnere m the last relation the RWA has been used. 
The almost degenerate next higher lying excited state |e, 0) is separated from 
the ground state by 002 = Ao + ^_^' 2 ~ A + A ®_ Q - These two frequencies 
show up in the spectrum, see Fig. 3. The perturbatively calculated values of 
the peak positions are for these parameters lo\ = 0.73A and cu 1 = 1.02A . 
The spectral weights of the peaks are different, i.e., cq 
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5.3 Positive detuning Q > A 

For the case of positive detuning il > A , also two frequencies appear, see 
Fig. 3. A perturbative treatment as in Subsection 5.2 gives the frequency of 
the peaks at ^ = A + « A + ^ and u 2 = Q - |f^_ ~ ft _ 

where in the last relations the RWA has been used again. For the parameters 
used in Fig. 3, the results are lo\ = 0.98 and uj 2 = 1.27 which agree well with 
the positions of the peaks in the spectrum. The spectral weights are in this 
case inverted, i.e., a± ^> a 2 - 

5.4 Dephasing rates 

In order to extract the dephasing rates Tj, we fit two Lorentzians to the peaks 
in the spectrum shown in Figs. 1 and 3. The dephasing rates follow according 
to Eq. (12). The dependence of the dephasing on the frequency of the HO is 
shown in Fig. 4. The main figure displays the result if the interaction strength 
g between TSS and HO is kept fix, while its inset depicts the result for a 
kept fix, i.e., g also varies with varying Q. Each figure displays the dephasing 
rates of both peaks. The dominating peak with the dominant spectral weight 
is peak 2 for negative detuning Q < A , while it is peak 1 for positive detuning 
Q > A . Both peaks show a clear maximum at which the dephasing is largest. 
This happens naturally in the vicinity of the resonance Q ^ Aq, where the 
effect of the coupling of the TSS to the HO is most effective and therefore 
dephasing of the HO penetrates through to the TSS best. 
The inset of Fig. 4 can be readily compared to the result of the recent work 
by Kleff et al, see Fig. 3b of Ref. [15]. There, the dephasing rate shows a 
particular almost discontinuous dependence on the ratio of fi/A around 1. 
In contrast to their findings, we do not observe such a type of behavior but 
instead a rather smooth dependence. The almost discontinuous feature is likely 
to be due to the fact that in Ref. [15] , the zero temperature dynamics has been 
considered. 



6 Varying the TSS-HO coupling strength g 

The dependence of the dephasing rates on the coupling g between the TSS 
and the HO is depicted in Fig. 5 for the case of negative detuning f2 < A . As 
expected, it increases monotonically with increasing g and for this case, peak 
2 is dominating, see the spectral weights in the inset of Fig. 5. When the 
coupling g between the TSS and the decoherent HO is increased, dephasing 
of the TSS is also larger. In general, the dependence of the dephasing on g is 
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Fig. 4. Main: Dephasing rates corresponding to peak 1 and peak 2 in the Figs. 1 
and 3 as a function of the HO frequency The parameters are: e = 0, g = 0.07Ao, 
k = 0.014, ksT = 0.1/iAo. Inset: Same for stronger damping k = 0.02 with 
a = 0.01 =const. (like in Ref. [15]). This implies that with varying Q also g is 
changed. 
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Fig. 5. Main: Dephasing rates in dependence of the coupling g between the TSS 
and the HO for negative detuning f2 = 0.75Ao- For the remaining parameters see 
Fig. 4. Inset: The corresponding spectral weights a«. One observes that peak 2 is 
dominating for this detuning. 



observed to be rather weak, as long as g < A , which is usually the case. 



7 The role of temperature T 

As in the pure Ohmic spin-boson system for weak damping, we also expect a 
rather weak dependence of the dephasing rate on temperature in the region of 
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Fig. 6. Temperature dependence of the dephasing rates Ti and T 2 for the same 
parameters as in Fig. 4 (main) for 17 = 0.75. Inset: Corresponding spectral weights 



low T for this more complicated system. Indeed, the results confirm this, see 
Fig. 6. The dephasing rates related to both peaks are almost constant with 
temperature T. 



8 Dynamics in presence of a finite bias e ^ 

In this section, we consider the case of a biased TSS. Due to this bias, the 
spectrum shows two more resonances: (i) the relaxation peak at to — 0, (ii) an 
additional dephasing peak. The frequency of the latter can be estimated by 
second order perturbation theory in g. It corresponds to the energy difference 
between the states |e, 0) and \g, 1), which is cj 3 = A b — il + . A typical 
spectrum is shown in Fig. 7. For this set of parameters, we obtain = 0.46A , 
which agrees well with the position of the numerically obtained resonance. The 
remaining two dephasing peaks at lo\^ correspond to those of the unbiased case 
with the substitution A — > A b . 

The relaxation and the dephasing rates, r r and T^, respectively, for increasing 
bias are shown in Fig. 8. Note that we show only the dephasing rate of the 
dominant peak 2, i.e., of that one with the largest spectral weight. The qual- 
itative behavior of the rates corresponding to the remaining two dephasing 
peaks is similar (not shown). As one observes, the decay rates decrease with 
increasing bias for this case. This can be understood by considering the posi- 
tion of the resonance at Q relative to the TSS level spacing A for zero bias 
in the spectral density J(w), Eq. (2). For the case shown in Fig. 8, Q < A . 
If we increase the bias e, the effective TSS level spacing A b increases. Since 
the spectral density is then decreasing (we move on the right flank of the 
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Fig. 7. Fourier transform S z (u) for the biased case e = 0.5. Compared to the 
unbiased case, two more resonances appear: the relaxation peak of width T r at 
<jj = and a third dephasing peak of width r3 at u = uj^. Parameters are: 
n = 0.75A , g = 0.13A , /c = 0.014(-> a = 0.003), /c B T = 0.1hA . 
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Fig. 8. Relaxation rate T r and dephasing rate T2 for varying bias e. Only the 
dominant dephasing peak is considered. The qualitative behavior of the other 
two rates is similar (not shown). Parameters are: f2 = 0.75Ao,g = 0.13Ao, 
k = 0.014(— >• a = 0.003), ksT = 0.1/jAq. Inset: Corresponding spectral weights 



a,-. 



resonance in Eq. (2) towards larger u), dephasing is less effective leading to 
decreasing decay rates. This reasoning can directly be verified by considering 
the opposite situation Q > A (not shown). In this case, an increasing bias 
moves the TSS level splitting towards increasing spectral weights, since we are 
moving on the left flank of the Q- resonance in Eq. (2), causing an increase of 
the dephasing and relaxation. 
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Fig. 9. Dephasing rates Ti and T2 in comparison with the weak-coupling rate 
for varying Q. Parameters are: e = 0,g = 0.07Ao, n = 0.014(— > a = 0.0005 for 
Q, = Aq), ksT = O.IHAq. Inset: Corresponding results for varying g for Q = 0.75Aq. 



9 Approximate analytical approaches 



The standard analytical approach to the reduced dynamics of a TSS weakly 
coupled to an environment consists in solving a master equation for the re- 
duced density operator [22,1]. Path-integral methods [1], as well as the Bloch- 
Redfield formalism can for instance be used [23]. The spin-boson system has 
been investigated in great detail for the case of weak-coupling to an Ohmic 
bath, having a continuous spectral density. The most general description dates 
back to Argyres and Kelly [23] where also strong time-dependent driving has 
been included. Recently, the equivalence between the the path-integral ap- 
proach and the Bloch-Redfield theory has been demonstrated in Ref . [24] . 

In this section we shall apply a weak coupling approach to our problem, and 
compare with the exact numerical results presented in the previous sections. 
In view of the possible application to more complex systems, it is in fact 
important to understand for this simple problem, in which regime the pertur- 
bative estimate differs from the exact solution. In addition, the weak-coupling 
approach is widely used to estimate the dephasing and relaxation rates in 
solid-state qubits [4,6,25]. 

In the second part of this section, we will show that a more appropriate ap- 
proach to the problem of decoherence due to a structured environment consists 
in enlarging the Hilbert space of the quantum system in order to include the 
harmonic oscillator. 
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9.1 Comparison with the spin-boson model for weak damping 

For the spin-boson system the dephasing and relaxation rates can be evaluated 
analytically to lowest order in the coupling to the bath [1,18]. For the case of 
a spectral density showing Ohmic behavior at low frequency, like Eq.(2), they 
read 

_T r ne 2 2ak B T 




where A b = ^JAq + e 2 . The results for the varying parameters Vt,g,T are 
shown in Figs. 9 and 10 for the symmetric case. The dependence on the bias e 
is depicted in the inset of Fig. 10. We show only the results for the dephasing 
rates, the relaxation rate T r behaves in a qualitatively similar way. 
We notice that the weak-coupling rates Eq. (13) give a good estimate of the 
dephasing rate only if f2 3> A 6 . This result can be understood considering 
the physics of the model in the form given in Eq. (3). If we regard the TSS 
as coupled to an environment represented by the damped HO, where g is 
the coupling constant, the rates in Eq. (13) follow from the usual Bloch- 
Markov treatment. This latter approach is valid when gr c <C 1, where r c 
represents the typical correlation times of the weakly damped oscillator [22]. 
In the temperature regime considered r c may be identified with the damping 
rate of the HO, 7 _1 [1]. Thus, the perturbative approach holds when g <C 7. 
With the parameters of Fig. 9, this condition is satisfied for ^ 3> 0.8. 

Of course the weak-coupling rates in Eq. (13) fail in describing the dephasing 
rates if the TSS is close to or at resonance with the HO. For VI ss A , the 
weak-coupling rates in Eq. (13) overestimate the dephasing because there are 
coherent exchange processes between the TSS and the HO mode which are 
not captured in a weak-coupling approach. It is interesting to notice that for 
Q < A the damped HO behaves as an environment with long correlation 
times, and the weak coupling theory fails because of memory effects. The 
results in Eq. (13) underestimates the dephasing rate in the regime Q < A , a 
feature which has already been noticed for environments responsible for 1/f 
noise [26]. 

In fact, a weak-coupling approach is appropriate as long as the correlation time 
of the bath is sufficiently short (as it is the case for an Ohmic bath with small 
a). However, the presence of the peak in the bath spectral density induces 
bath correlations which decay only on a time scale given by the width of the 
resonance. 
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Fig. 10. Same as Fig. 9 for varying temperature T. Parameters are: 
£ = : Q = 0.75A ,c/ = 0.13A , k = 0.014(^ a = 0.003). Inset: Corresponding 
results for varying the bias e for k^T = 0.1^A . 

9.2 Three-level system as an approximative description 



In order to obtain an analytical expression for the dephasing rates, we restrict 
for simplicity to the simplest case of resonance Q = Aq between the TSS and 
the HO. Since the deviation from the pure Ohmic case is most pronounced 
here, this is also the most interesting case for an analytical solution. We fur- 
thermore assume an unbiased TSS, 5 = 0. All the following considerations can 
be generalized to the off-resonant case and also with a biased TSS. The ex- 
pressions, however, will become more involved. As it turns out, the TSS+HO 
can be restricted at low temperatures to a three-level-system (3LS) having 
performed second order perturbation theory in g. Assuming weak coupling 
to an Ohmic bath, we can derive the Redfield equations which can be solved 
analytically. By this, we obtain closed expressions for the dephasing rates Tj 
and the frequencies Ui. 



9.2.1 Redfield equations 

To start, we diagonalize the TSS+HO Hamiltonian which corresponds to the 
substitution o~ x — > o z and o~ z — > —a x (the resulting transformed operators are 
denoted by the overline). Moreover, we apply the rotating wave approximation 
(RWA) which is appropriate at resonance. Together with a x = a + + er_, we 
obtain H^^ +HO = — ^^-o- z + KIB^B — hg(Ba + + B^a). This is the Jaynes- 
Cumming Hamiltonian which can be diagonalized exactly. The groundstate 
|g, 0) is not shifted while the degeneracy between the eigenstates of the un- 
coupled system is lifted. The higher lying eigenstates appear in pairs where 
the energy of the n-th pair is nhVL above the ground state. Each pair is split 
by HA n = Ihg^pri. As already used in Section 5.1, we can restrict at low tem- 
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peratures and weak damping the system only to the ground-state \g, 0) and 
the first pair of excited states, which are denoted as 

|-, 0) = -L[| e , 0) + \g, 1>] , |+, 0) = ^=[-\e, 0) + \g, 1>] , (14) 



in the dressed states basis. The corresponding energies of this 3LS are E_ = 
H(Q — g) and £+,0 = h(Q + g). For our observable being the population differ- 
ence, this implies that P + (t) = P(t) = (a z ) t = tr TS s(p(t)cr z ) = tr TS s(p(t)cr x ). 
In the basis of the dressed states, we find 

Pit) = V2Re((g,0\p(t)\-,0) - (g, 0\p(t)\+, 0» . (15) 



Hence, we need to solve the Redfield equations for the off-diagonal elements 
P ff o-o(0 = (fi'>0|p(t)|-,0) and p fl0 ,+o(0 = 0|p(t)|+, 0). A standard deriva- 
tion of the Redfield equations [1] gives (to simplify the notation, we omit in 
the following the overline to denote the dressed states basis) 



p g0 ,+o = + g)p g o,+o - ^{C(n + g) + C(Q - g) + C*(-Q - g)}p g0 , +0 
~C*(-Q + g)p g0 ,-o + \{C(-n -g)+ C*(n + g)}p +0 , 9 o 
+ ^{C(-n + g) + C*(n + g)}p^ g0 

Pg0 ,-o = i(n - g)p g0 ,- - ^{C(n + g) + C(Q - g) + C*(-Q + g)}p g0 ,-o 
~CT(-n - g) Pg0 , +0 + \{C(-n -g) + C*(Q- g)} P+0 , g o 
+\{C{-n + g ) + C*(n- g)}p- , g o ■ (16) 
Here, C (u) is the one-side Fourier transform of the bath autocorrelation func- 



tion C(t) = / °° <L/ Johm («/) 



coth ^4r- cos uj't — i sin u't 



2 



. It follows that 



ReC(cu) = ^ J hm(^)[coth - 1] 



oo 



lmC(cu) = J dJ coth ^ - u'] = A(u) - A'(uj) . (17) 

o 

In a next step, we use C(±f2 ± g) m C(±Q) and perform the secular ap- 
proximation by neglecting terms in p^ \ for which the Redfield tensor elements 
\Rnm,ki\ ^ l^nm — ^ki\- In particular, this implies to neglect the terms with 
P+o, g o an d P-o, g o- We arrive at the final form of the Redfield equations 
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7i T2 

A?0,+0 = + g)PgO,+0 - yP<yO,+0 ~ yP fl O,-0 , 

7i T2 

p 9 o,-o = - #)p 9 o,-o - yP<?o,-o - yP 9 o,+o , (18) 

where the rate constants are given by 71 = 2C(Vt)+C*(— Q) and 72 = C*(—Q). 
This system of coupled first-order differential equations can be solved by di- 
agonalizing the coefficient matrix. The eigenvalues are given by 



A± = -y \/72 2 -V. (19) 



The dephasing rates follow as T± = Re\± and the oscillation frequencies are 
given by uj± = ImA±. 

The Redfield equations are expected to be applicable for weak damping, i.e., 
k C 1. Moreover, this 3LS approximation can only be applied as long as 
the temperature is low enough, i.e., /c#T <C ft£L. In addition, it is required 
that kMI <C ksT. To verify our approximation, we compare in the following 
subsections the results with the numerically obtained values of QUAPI. 



9.2.2 Dephasing rates 

Plugging in the Ohmic spectral density (4) in Eq. (17), the dephasing rates 
r± follow as the real parts of A± in Eq. (19). We find 



± 



3 coth - — — 1 



1 A ^ 

± -A cos - , 

2 2 



(20) 



where 



(p = arctan 





" co th tf. + 1 


I 




r cot h^ + l" 


2 - P - Ag 2 



(21) 



and 



A 



7T 



coth — - + 1 
2 



coth — - + 1 
2 



-i 2 - v 

1/4 



(22) 



Moreover, we have / = A(Q) + A'(fi), where 
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2nk B T 



+- 



( ^Eif-Ue^Eif--) 



(23) 



and 



e n/^ Ei ( 



(24) 



Here, ^(z) is the digamma function and Ei(z) denotes the exponential inte- 
gral. Choosing the following dimensionless parameter set: k = 0.0014, T = 
0.01, uj c = 10, g = 0.07H « = 0.00005), we obtain T+ = 0.0025 and T_ = 
0.002, which have to be compared with the QUAPI results I\ = 0.0030 and 
r 2 = 0.0027. The 3LS results deviate from the QUAPI results by 17% and 
25%, respectively. 



9.2.3 Oscillation frequencies 

We find for the oscillation frequencies 

u ± = n - -A(Q) + -A' (SI) ± -A sin ^ . 

2 2 2 2 



(25) 



For the same parameter set as above, we obtain uo + = 0.93 and c<j_ = 1.07 from 
the 3LS and from QUAPI, we find uj\ = 0.93 and 002 = 1.07. The difference 
between the two rates is Au = 2g, as accurately predicted by the result from 
the exact diagonalization of the Hamiltonian in RWA. 



10 Application to experimentally realized superconducting flux 
qubits 



In this section, we apply our model to qubit devices which have been realized 
experimentally in the Delft group of J. Mooij. We refer to two devices, i.e., (i) 
the more recent flux qubit of Ref. [6], where the SQUID was directly attached 
to the qubit in order to increase the mutual coupling, and (ii) the flux qubit 
of Ref. [4] which was inductively coupled to a surrounding dc-SQUID being 
not directly in contact. 

In order to extract the relevant parameters for our model, we use the cor- 
respondence of the qubit-SQUID setup with an electronic circuit [25]. This 
allows us to express the parameters in our model in terms of the experimen- 
tally available parameters. 



20 



Fig. 11. Sketch of the circuit model of the dc-SQUID consisting of an ideal induc- 
tance Lj which is shunted by the capacitance C s . The Ohmic resistance is captured 
by Ri. For more details, see Ref. [25]. 



In this approach, the damped dc-SQUID provides an electromagnetic envi- 
ronment for the qubit. It can be characterized by the impedance of the corre- 
sponding circuit sketched in Fig. 11. 

The dc-SQUID far away from the switching point is modelled by an ideal 
inductance Lj. In the investigated devices, a large shunt capacitance C s was 
designed in order to control the environment. All Ohmic resitance of the circuit 
is modelled by R[. Together with the flux quantum the plasma frequency 



of the SQUID can be calculated as uj p = i^ jC = y 2n f^j~ ^ • Here, I c q is 
the critical current of a single Josephson junction (the SQUID possesses two 
junctions) and / = 7r$ e xt /^o is the scaled external flux, which is also denoted 
as the frustration of the SQUID. By comparing our expression (2) for the 
effective spectral density with Eq. (7) in Ref. [25], we find for our parameters: 



1 1 {2eI p I b M q f 1 

n = W " K= CJt^ a= C^ Rl 2< (26) 



For the coupling between qubit and SQUID, the persistent current I p of the 
qubit interacts via the mutual inductance M q with the bias current i& of the 
SQUID. 

An alternative way to obtain the Hamiltonian (3) has been worked out by H. 
Nakano and H. Takayanagi [27] by starting from a microscopic Hamiltonian. 
It carefully takes into account the phases for the three Josephson junctions in 
the qubit ring and also the two junctions of the SQUID. By integrating out 
the corresponding hidden degree of freedom and by approximating the SQUID 
potential well for low bias currents by a harmonic oscillator, one arrives at a 
Hamiltonian equivalent to our Eq. (3). 
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10.1 Qubit 1 



For the Delft qubit of Ref. [6], the measured value of the level spacing is 
A = 2ir 3.4 GHz. The persistent current of the qubit is I p = 330 nA. The mea- 
surements of the dephasing and relaxation times were done for the biased qubit 
with an effective level separation A& = 27r5.71 GHz. The dc-SQUID is spec- 
ified by the parameters I c0 = 2.2fiA,C s = 10 pF and / = arccos(|) = 1.45. 
The plasma frequency follows as lu p = 27r2.0 GHz. 

In the Delft experiment, the inductive coupling between the qubit and the 
SQUID was rather large, i.e., M q = 9 pH and the circuitry environment is 
assumed to have an Ohmic resistance Ri ~ 160f2. The measurement time is 
therefore short. Hence, the experiment detected dephasing and relaxation of 
the qubit during operation during which the qubit is in principle decoupled 
from the detector. However, due to the experimental design, the decoupling 
is not perfect. A rather likely small asymmetry between the two Josephson 
junctions in the SQUID generates a small circulating current which is com- 
pensated by a small bias current flowing through the SQUID. Assuming [28] 
that this bias current is about 5% of the critical current J c0 of a junction, this 
yields an effective bias current also during operation time of I b — 110 nA. The 
experiment has been performed at a temperature of T = 25 mK. 
Plugging in the parameters and scaling them with respect to the pure level 
spacing Ao, we arrive at the dimensionless parameters e = 1.349, Q = 0.59, k = 
0.0078,# = 0.008, T = 0.15, (a = 1.1 x 10~ 5 ). We assume for the cut-off fre- 
quency uj c = 10. With similar QUAPI parameters as above (At = 0.06, M = 
12 and K = 2), we find for the dephasing rate = 0.0001A , which implies 
a dephasing time of = 470 ns. It deviates by more than a factor of 20 from 
the measured value of T^exp = 20 ns. 

In order to check the applicability of the weak coupling approach, we cal- 
culate the dephasing rate for the above parameter set. The weak-coupling 
dephasing time is calculated to be r d iWC = 6.6 /is. The weak coupling approach 
underestimates dephasing in this case by more than one order of magnitude. 
The determination of the relaxation rate for this specific parameter set is not 
possible with QUAPI since we would have to iterate the dynamics up to un- 
realistically long times. 

The relatively large calculated dephasing and relaxation times present good 
news for a possible future use of the qubit in a quantum computer. We have 
assumed as noise source the Johnson-Nyquist noise from the electromagnetic 
environment. The reason for this rather large deviation from the experimen- 
tal data is that in the experiment an accidental additional resonance with an 
environmental mode occurred. This mode stems from a superconducting loop 
designed around the full qubit-plus-SQUID device. Its resonance frequency 
has been determined experimentally to be around 3 GHz, which is close to 
the qubit level spacing at the degeneracy point. This kind of resonance will be 
avoided in the design of the next generation of the qubit which is currently on- 
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going [28]. An additional filtering of external electronic noise amplified by the 
various measuring devices is also expected to increase the measured dephasing 
time. 

10.2 Qubit 2 

We have also calculated the dephasing rate due to Johnson-Nyquist noise 
for the flux qubit reported in Ref. [4]. The parameters for this qubit are 
A = 27r0.66 GHz and I p = 500 nA. We determine the dephasing rate at 
the degeneracy point, i.e., e — 0. The dc-SQUID is specified by the parame- 
ters 2J c0 = 200nA, C s = 30 pF and / = 0.75. The plasma frequency follows as 
u p = 2tt0.61 GHz. 

The device was designed in the way that the dc-SQUID enclosed the qubit. 
The mutual inductance is M q = 8 pH and the Ohmic resistance of the leads 
is estimated as Ri = 10011 In this design, the dc-SQUID dominates the de- 
phasing of the qubit since the mutual coupling is weak and the measurement 
required to average over many switching events of the qubit. The bias current 
of the SQUID is taken to be close to the switching current, i.e., h = 120 nA. 
The experiment was performed at T = 30mK. 

The dimensionless parameters follow again from scaling with respect to the 
level spacing A . We find SI = 0.93, re = 0.014, g = 0.063 (a = 5 x 10" 4 ) and 
T = 0.95. The dephasing rate is determined by QUAPI with the same param- 
eters as above for qubit 1. We find Yd = 0.05A yielding a dephasing time of 
Td = 5ns. Also in this case, we can compare this result with the outcome of 
the standard weak-coupling approach which yields the value of = 5 ns. The 
measured value for the dephasing time is Td >exp = 10 ns. 



11 Conclusions 

In summary, we have investigated the dynamics of the spin-boson problem for 
the case when the frequency distribution of the bath shows a distinct reso- 
nance at a characteristic frequency Q. We have mapped this model onto that of 
a two-state-system (TSS) coupled to a single harmonic oscillator (HO) mode 
with frequency Q, the latter being weakly coupled to an Ohmic bath. Since 
an Ohmic bath induces rather fast decaying correlations at finite tempera- 
tures, the numerical method of the quasiadiabatic propagator path-integral 
(QUAPI) can be applied by using the TSS+HO as the central quantum sys- 
tem exposed to Ohmic damping. This becomes treatable by QUAPI as the 
later relies on the technique of cutting the correlations after a finite memory 
time. Having then studied the decay rates for dephasing and relaxation of a 
quantum superposition of eigenstates of the TSS, we found new features of the 
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dynamics. This includes at least two inherent frequencies of the time evolution 
of the TSS, as well as pronounced increased dephasing and relaxation if the 
TSS is in resonance with the localized harmonic mode. Our findings clearly 
demonstrate that the whole frequency spectrum of the peaked environment, 
and not only its low frequency Ohmic part may become relevant in deter- 
mining the TSS relaxation dynamics, if the frequencies of the TSS and the 
HO are comparable. In particular, we showed that, for a peaked environment 
close to resonance conditions, the commonly used estimate of the decay and 
relaxation rates to lowest order in the coupling between the TSS and the har- 
monic reservoir becomes inadequate. The appropriate route is to consider the 
TSS+HO as the quantum system, and then perform perturbation theory in 
the coupling between such a system and the smooth Ohmic environment. To 
this extent, we established the Redfield equations for the unbiased, resonant 
TSS-HO system. It turns out that this system can be reduced to a three-level 
system at temperatures well below the characteristic frequency Q. Within 
the three-level system approximation, we have obtained analytic expressions 
for the relaxation and dephasing rates, which well agree with the numerical 
QUAPI rates. 

Finally, we have applied the general model to experimentally realized super- 
conducting flux qubits coupled to a dc-SQUID playing the role of the detector. 
We find considerably long dephasing times which represents an encouraging 
result for a possible use of these devices as qubits for future quantum infor- 
mation processors. 
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